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Abstract: Optical flow is the apparent (or perceived) motion of image brightness pat- 
terns arising from relative motion of objects and observer. Estimation of the optical flow 
requires the application of two kinds of constraint: the flow field smoothness constraint 
and the brightness constancy constraint. The brightness constancy constraint permits 
one to match image brightness values across images, but is very restrictive. We propose 
replacing this constraint with a more general constraint, which permits a linear transfor- 
mation between image brightness values. The transformation parameters are allowed to 
vary smoothly, so that inexact matching is allowed. We describe the implementation on 
a highly parallel computer, and present sample results. 

1 Introduction 

Optical flow is the apparent (or perceived) motion of image brightness patterns arising 
from relative motion of objects and observer. Optical flow can give important information 
about motion of the observer (i.e. passive navigation), motion of objects in the scene, and 
the spatial arrangement of these objects. Additionally, discontinuities in the optical flow 
field can be used to segment the image into regions corresponding to different objects. 

To be precise, an optical flow field is a two-dimensional vector field relating brightness 
patterns in an image at one instant of time to brightness patterns at the next instant of 
time. There does not exist a unique optical flow field for a given image sequence; rather 
there are infinitely many flow fields satisfying the image constraints. This illustrates one 
of the difficulties associated with determining optical flow, namely, identifying sufficient 
constraint to produce a unique optical flow field. The other difficulty in determining 
optical flow is more fundamental, and involves finding image elements to be placed into 
correspondence. 

As in stereo, methods for computing optical flow can be classified according to whether 
detected features are used as primitive elements, or whether image brightness values (and 
gradients) are used directly. Feature-based approaches to optical flow use detected edges 
almost exclusively (Hildreth [1983], Davis et al [1983], Murray & Buxton [1984], Wohn 
[1984]). Although it is possible to use detected points to process visual motion, no one 
seems to have attempted to determine optical flow from isolated point displacements. 
Horn & Schunck [1981] were among the first to use image brightness directly to determine 
optical flow. They solved the problems mentioned above by identifying two constraints: 
the spatial smoothness constraint and the brightness constancy constraint. Other methods 
based on various smoothness assumptions have been proposed (Prager & Arbib [1983], 
Paquin & Dubois [1983], Yashida [1983], Anandan [1984]). 

The smoothness constraint arises from the observation that most visual motion is 
the result of objects of finite size undergoing rigid motion or deformation. Neighboring 
object points have similar motions or velocities, and to the extent that they project 
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to neighboring image points, neighboring image points will also have similar motions. 
Therefore, the optical flow field should be smooth almost everywhere. Exceptions occur 
at occluding boundaries, where neighboring image points are not generally the projections 
of neighboring object points. 

The brightness constancy constraint rests on the assumption that the brightness of 
a small image patch remains approximately constant as the corresponding surface patch 
moves in the environment. This is a reasonable assumption when the lighting conditions 
are unchanged between successive images, object surfaces are non-specular, and there is 
only a small amount of motion between image frames. If these conditions are met, then 
the brightness constancy constraint will apply approximately at all image points. 

It should be noted that the smoothness constraint depends on the scene structure, 
and is independent of illumination, surface reflectance characteristics, and the type and 
degree of motion involved. On the other hand, the brightness constancy constraint does 
not depend on scene structure (except for the influence of surface microstructure on 
reflectance), but instead depends on the degree and types of motion, and factors such 
as illumination and surface reflectance which affect image irradiance. Violations of the 
brightness constancy constraint, when they occur, affect image patches or even entire 
images. Therefore, it is important to find ways to relax this constraint. 

Cornelius & Kanade [1983] propose a variation of the Horn & Schunck [1981] method. 
In their formulation, they allow gradual changes in the way an object appears in a se- 
quence of images. An image point does not have to preserve the same brightness value 
as the object point that give rise to it moves in the environment, however, the variation 
is enforced to be smooth from one image point to the next. 

In this paper, we propose a new formulation by relaxing the brightness constancy con- 
straint. Our approach does not require exact brightness matching across image frames, 
but accepts even approximate matches. We achieve this by permitting a linear transfor- 
mation of image brightness values between image frames, and constraining the allowed 
transformations. Our formulation, in special cases, reduces to that of Cornelius & Kanade 
[1983] or that of Horn & Schunck [1981]. 

2 Mathematics of Optical Flow 

Let a coordinate system be aligned with the imaging system so that the 2-axis points 
along the optical axis. The image plane can arbitrarily be chosen to lie at z = 1 so that 
image points are given by r = (x,y, 1) T . Let E(r, i) denote the brightness of image point 
r at time t. At a later time t + 6t, the brightness pattern at r will have moved to a new 
location r + 6r = (x + 6x, y + 6y, 1) T . The optical flow is the velocity field arising from 
the perceived motion of image points r. It is derived from the displacements of image 



points by taking the limit as 6t — » 0, 

dr ,dx dy . T . .</• 

where u and v denote the components of optical flow r* . 

The brightness constancy constraint of Horn & Schunck [1981] expresses the restriction 
that the brightness of an image patch remains approximately constant as the surface patch 
that gives rise to that image patch moves in the environment. Setting the total derivative 
of image brightness equal to zero, we can write 

dE 



, =0- 
dt 



Applying the chain rule, we obtain 



dE dEdx dEdy n 

1 1 = 0, 

dt dx dt dy dt 

or 

E t + E r • r t = 0, 

where E r = (dE/dx,dE/dy,0) T . This equation is sometimes referred to as the image 
brightness change constraint equation derived under the constant brightness assumption. 
It has also been referred to as the optical flow constraint equation. 

Assuming brightness constancy, we define optical flow as any 2D vector field r*, defined 
on the image plane, that satisfies the image brightness continuity equation. 

2.1 Image Brightness Constraint 

Nagel [I983a,1983b] suggests a formulation that incorporates second-order effects in or- 
der to obtain a better estimate of the optical flow around edges and corners. In his 
formulation, the brightness change constraint equation is written 

E t + E t • r t + -rjE rr r t = 0, 

where „ 

(d 2 E/dx 2 d 2 E/dxdy 0' 
d 2 E/dxdy d 2 E/dy 2 
0. 

Here, again, the constraint equation rests on the brightness constancy assumption. 

The formulation proposed by Cornelius &; Kanade [1983] allows gradual changes in 
the way an object appears in a sequence of images. In their formulation, the brightness 

change constraint is written 

dE 

— = E t + E r - r t . 

dt 
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An image point does not have to preserve the same brightness value as the object point 
that gives rise to it moves in the environment. Hence, the rate of brightness change can 
be non-zero; that is, 

dE ^ n 

More generally, we propose a formulation that allows a linear transformation between 
brightness values in consecutive images. We choose a linear transformation because it is 
one of the simplest non-trivial transformations. This is a less restrictive assumption than 
brightness constancy, and can be formulated as 

E(r + 6r,t + St) = M(r, t) E(r, t) + C(r, t) , 

where M is the multiplier and C is the offset functions in the linear transformation. This 
is our revised image brightness change constraint equation. 

For small St, we expect M to be close to 1, and C to be close to 0. Since we are 
dealing with incremental changes, we can let M = 1 + Sm and C = Sc. In fact, m and c 
are the quantities of interest to us. Noting that Sm — *• and Sc — ► as St — *■ 0, we can 
define time derivatives, mt and ct, 

Sm Sc 

m* = lim -t— and c* = lim — , 

6t->0 St 6t->0 St 

that we will use in our derivation. 

Rewriting the brightness change constraint equation, we obtain 

E{r + Sr,t + St) = [1 + Sm(r,t)} E(r,t) + Sc{r,t). 

The left hand side can be expanded as follows: 

dE dE 

E(r + Sr,t + St) = E{r, t) + — -6r+ —St + 0{e) = E r -Sr + E t St + O(e). 

or at 

Substituting this into the constraint equation and simplifying, we have 

E r • Sr + E t St - E Sm - Sc + O(e) = 0. 

Finally, dividing through by St and taking the limit as St —*■ 0, we arrive at 

E t + E, • r t - E m t - c t = 0. 

This is our revised optical flow constraint equation. In the special case that M = 1 and 
C/0 (and, hence, mt = but ct ^ 0), this becomes similar to the constraint equation 
in the Cornelius & Kanade [1983] formulation: 

ct = Et + E r • Tf. 

Further, in the more restricted case that M = 1 and C = (that is, mt = Ct = 0), our 
constraint equation reduces to the one in the Horn & Schunck [1981] formulation: 

E t + E r • r t = 0. 





Figure 1. Corresponding iso-brightness contours in an image sequence. 

3 The Aperture Problem 



There are an infinite number of valid optical flows given an image sequence. To see 
this, we note that if Tt is an optical flow, so is 

r' t = r t + f{x,y) {E r x z) 

for any f(x, y), where z is a unit vector perpendicular to the image plane. First, note 
that r' t • z = rt • z = 0, as it should. Furthermore, we have 

E t + E T - [T t + f {E r xz)}-Em t -c t = E t + E r -rt-Em t -ct = 0. 

This can also be explained graphically. Consider the simple case where M = 1 and C = 0, 
so that the constraint equation reduces to that of Horn & Schunck formulation: 

E r • r t + E t = 0. 

Referring to figure 1, suppose that C" is a contour of constant brightness in the second 
image corresponding to contour C in the first image. It is not easy to decide which point 
P' on C corresponds to a particular point P on C since the contour generally changes 
shape as the object moves in the environment (Horn [1986]). In fact, there are many 
possible ways to establish correspondence between points on contours C and C' . This 
ambiguity has been referred to as the aperture problem. In terms of the iso-brightness 
contours, any vector field that transforms contour C into contour C is an acceptable 
optical flow. 
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In our extended formulation, the indeterminacy of the optical flow is even worse than 
this example suggests. Since m< and ct are unconstrained, the optical flow field can be 
completely arbitrary, with either or both of these transformation fields varying in such a 
way as to guarantee that the brightness constraint is obeyed. It is, therefore, necessary 
to select, out of the infinite number of possible optical flows, one which is consistent with 
the physical constraints of the problem. One may hope to obtain an optical flow field 
that approximates the apparent motion of brightness patterns in the image as judged by 
a human observer. 

3.1 Smoothness Assumption 

Discontinuities in depth (for example, at occluding boundaries) give rise to disconti- 
nuities in the optical flow field. Also, object motions may be different across occluding 
boundaries, which can give rise to discontinuities in the optical flow. Additionally, discon- 
tinuities can be expected in mt and ct, if illumination conditions or reflectance properties 
that depend on surface material change abruptly as the surface moves in the environment. 
In the absence of depth discontinuities or abrupt changes in illumination or surface re- 
flectance properties, the optical flow and transformation fields are expected to be smooth. 
Based on these facts, we require that the optical flow, the multiplier, and the off-set fields 
should be consistent with our revised optical flow constraint equation, and should vary 
smoothly from one image point to the next. 

Smoothness can be imposed by minimizing a functional that is a measure of departure 
from smoothness. Horn & Schunck [1981] proposed minimizing the integral of the square 
of the magnitude of the gradient of the optical flow. Hildreth [1983] investigated a similar 
formulation, but incorporated different measures of smoothness. 

The gradient of the optical flow is 

~ / du/dx du/dy 

Vr* = -7— = ( dv/dx dv/dy 

dr \ 



The measure of departure from smoothness that is to be minimized is written 

12 dxdy. 



= //|Vr,|| 



Here, || • \\* denotes the Euclidean or Frobenius norm of a matrix, which is the sum of the 
square of all the elements of the matrix. 

Similarly, smoothness deviations can be defined for the transformation fields 

e m= \\Vrnt\\l<lxdy and e c = I I \\Vc t \\\dxdy. 



4 Minimization 

The image brightness constraint and the smoothness constraints can be combined by 
defining a single functional that weighs each contribution. Rather than enforcing the 
brightness change constraint exactly, we use a penalty term that measures the square of 
the error in the constraint equation over the whole image: 



«,-//(* 



t + E, - rt — E mt — c*) 2 dx dy. 

To ensure that the optical flow and the transformation fields (approximately) satisfy the 
optical flow constraint equation, we want e& to be small. 

All together, the problem can be formulated as that of minimizing the functional 

C = ej -|- ng6g -f- ^m^m T AfiCc, 

where \ 8 , X m , X c weigh the total error contributed by each term. 

Using variational calculus, the Euler-Lagrange equations for this problem can be 
found. These equations form a set of necessary conditions that a solution to our mini- 
mization problem has to satisfy. Sufficiency is not guaranteed, in particular, it is possible 
for a particular proposed solution to obey the Euler-Lagrange equations yet not be a 
global minimum. This will occur at local minima, points of inflection, and local maxima. 
Note, however, that there is no global maximum, as the functional is not bounded from 
above. We will not address the question of sufficiency further in this paper. 

The variational problem is solved by using the formula 

where \£ is the integrand in the cost functional and / is each of u, v, mt or ct, in turn. 
Applying the above formula, we obtain 

E 
V 2 u = -^-{E t + E x u + E y v-Emt- c t ), 

E 
V 2 « = -^-{E t + E x u + E y v-Em t - c t ), 

p 

V 2 m t = -T—{E t + E x u + E y v-Em t - c t ), 

Am 

V 2 c t = ^-(E t + E x u + E y v-Em t - c t ). 

For a well-posed problem, we need to specify the appropriate boundary conditions. In the 
absence of fixed boundary conditions (values of u, v, mt, and ct on the image boundaries), 
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we need to specify natural boundary conditions. For our problem, the natural boundary 
condition is 

(/*,/s,) T -fi = 0, 

where n is a unit vector perpendicular to the boundary. Again, / can be any one of 
u,v,mt or Cf, in turn. 

4.1 A Discrete Implementation 

In the discrete domain, the Laplacian operator V 2 can be approximated as a center- 
surround operator 

v 2 /««(7-/), 

where / (the "surround") is an averaged or smoothed version of /. For example, we can 
use the following simple approximation: 

/ = j(/t'+l,y + fi-l,j + fij+l + fi,j-l)- 

The scale factor k can be absorbed into the appropriate A and, therefore, need not be 
considered further. 

Substituting the approximation to the Laplacian in the Euler-Lagrange equation, 
derived earlier, we can write a single matrix equation 

Af = g(f), 



where 



and 



f = 




E X Ey 



g(f) = 



X 8 u-E x E t \ 
\ s v-E y E t 
X m rh~i + EEt 
\ c ci+E t ) 



[E x + \ g E x E y 



E 2 y + X, 



-E X E 
■E y E 



-E X E -E V E E 2 + \ r 



-E x \ 

-Ey 

E 



V -E x 



-E v 



E 



;y ~ 1+XcJ 

These equations have to be solved iteratively since the optical flow, at each image cell, 
depends on the average of the optical flow from the neighboring cells. The same is true 
for mt and c*. 

Solving for the unknown fields, «, v, mt, and ct, we find that 



f = A- 1 g(f), 



where 



a 



and 



A c A m A,+A ro A,+ TP TP \ \ TP TP\ \ TP \ \ 

TP TP \ \ A,,A m A,+A m A,+ TP TP\ \ TP \ \ 

-rj x £jyA c A m E*XcX.+E*X c X m £jy^A c A 6 £j y A m A a 

E x E\ c \ a EyEX c X 8 (El+E>'+x.)x c x. ~E\\ 

E x \ m X s E v X m X e -EX S (E*+E*+'x.)x m x, J 



oc = X m X\ + E 2 X c X 2 g + {El + E 2 y + X s )X c Xm*s 



This is a system of linear equations, which can be used to recover the optical flow u and 
v, and the transformation fields mt and c<. 

The field equations can be solved iteratively, at every image cell, according to the 
equation 

f^A-'gCf*), 

where k is the iteration number. The matrix A (or A -1 ) depends only on the observed 
data. It needs to be computed once, but it differs from point to point. 

This formulation, in general, requires a lot of computation and is not really suitable 
for implementation on a serial machine. It can, however, be readily implemented on a 
highly parallel computer, such as the Connection Machine™. For a 128 by 128 image, the 
Connection Machine implementation runs approximately 1000 faster than a Symbolics 
3640™ Lisp Machine implementation. 

5 Examples 

Example 1 - Multiplier Effect with No Offset: Figure 2 shows a pair of image 
frames from a synthetic motion sequence. Each image contains a background texture 
and a central texture; each texture is gaussian-smoothed uniform noise. Sharp texture 
boundaries between the background and the central circular object have been preserved 
in each image. The circular region undergoes rigid rotation between frames. In addition 
to the rotation of the central circular region, the second image has been multiplied by 
a factor which varies linearly from 0.75 in the lower left corner to 1.25 in the upper 
right corner, as shown in Figure 3a. An offset of zero was used and all A parameters 
were set to 1. The computed multiplier field, shown in Figure 3b, varies in a range from 
0.82 in the lower left corner to 1.22 in the upper right corner. Here the linear trend is 
clearly discernible. The computed optical flow after 100 iterations is shown in Figure 
4 (the needles indicate flow direction and rate). The offset field for this example had 
negligible values, with absolute values not exceeding 0.002. This is to be expected, as the 
experiment was designed so that the offset field would not be needed. 
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Figure 2. A pair of images from a synthetic motion sequence; the central circular region undergoes rigid 
rotation. In addition to the rotation of the circular region, the second image has been scaled by a factor 
that varies linearly from 0.75 in the lower left corner to 1.25 in the upper right corner. 
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Figure 3. The true (a) and the computed multiplier fields used for the image sequence in Figure 2. 
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Figure 4. Example 1: The computed optical flow using the method described in this paper (all A param- 
eters were set equal to 1). 
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Figure 5. Example 1: The computed optical flow using the method described in this paper (A p = 0.1, 
but the remaining A parameters were set equal to 1). 
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Figure 6. Example 1: The computed optical flow using the method of Horn k Schunck. 
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Figure 7. Example 1: The computed optical flow using the method of Cornelius k Kanade. 
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More accurate results were obtained for X p = 0.1 (all other A parameters were set to 1) 
as shown in figures 5. The off-set field was negligible with absolute values not exceeding 
0.0002, and the multiplier field varied in a range from 0.76 in the lower left corner to 1.26 
in the upper right corner. 

To check the improvement offered by this algorithm, the Horn & Schunck algorithm 
was also used on this image sequence. This is equivalent to using A m = A c = oo in our 
formulation. The results of the Horn &; Schunck algorithm, Figure 6, were in agreement 
near the center of the image and the upper left and lower right corners, where the multi- 
plier was approximately 1. The two algorithms did not agree, and the unmodified Horn 
&Schunck algorithm was clearly incorrect, at the lower left and upper right image corners 
where the multiplier had its greatest effect. This illustrates the inability of the Horn & 
Schunck algorithm to correctly handle images sequences where the brightness constancy 
constraint does not apply. 

Figure 7 shows the solution obtained using the algorithm of Cornelius & Kanade. 
(This was done using A m = oo and A c = 1 in our formulation.) As expected, there is 
not much improvement over the solution from Horn & Schunck algorithm since their 
algorithm is designed to compensate for effects similar to an offset in an image sequence 
(however, the offset was set to zero for this example). 

Example 2 - Multiplier and Offset Effects: Figure 8 shows the pair of images for 
this example. The motion is as in the previous case, the multiplier field varies linearly 
from 0.9 in the lower left corner to 1.1 in the upper right corner, and the offset is 5 (the 
grey-level values were increased by 5 units). The computed optical flow using our method 
is shown in Figure 9 (all A parameters were set equal to 1). Figures 10 and 11 show the 
same using the methods of Horn & Schunck and Cornelius & Kanade, respectively. Again, 
these were obtained by setting A m = A c = oo (to simulate Horn & Schunck algorithm) , 
and A m = oo and A c = 1 (to simulate Cornelius & Kanade algorithm) in our formulation. 
These results are reasonable where the multiplier and offset effects approximately cancel 
each other. They break down where the multiplier and offset have their greatest effects; 
that is, in the upper right corner of the image. We conclude that the Horn & Schunck 
algorithm does not correctly handle images sequences where the brightness constancy 
constraint is violated. Similarly, the algorithm of Cornelius & Kanade breaks down in 
regions of the image where multiplier effects are dominant. 

6 Reducing the Computation 

In many practical situations, the off-set a is near 0. Therefore, the algorithm can be 
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Figure 8. A pair of images from a synthetic motion sequence; the central circular region undergoes rigid 
rotation, the multiplier field varies linearly from 0.9 in the lower left corner to 1.1 the upper right corner, 
and the offset is 5 units. 
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Figure 9. Example 2: The computed optical flow using the method described in this paper. 



19 



I 1 .. 



4 ,_^>.=^^- 



i :: J: 
± ;:f::<: 



:'"::/; 



-\ V 



i:./ 






/ x ? ; : /V 



:/:-.::/: 






:t.: :: :l . : :v 



■--./: 



:./-. : : i 



^H ; ;S^ : :*:^i^n 



--:: 7 ,:, ; ;;,::-::::^:::::: 
i i/::^: :/:::/!: ::|:::i: :/::,"::: :!:: 

:::::: \ : : v :::;?;::::.:■■■::::::::::: : 
■.., V ■< \ -, ~-- '*- i --^ 

Elii 



:::■::[:: 



■ - 1 > v ■% ■ 

::::::::: :V: : :\: :: :V : :t : 



: :::\: 



■!::V: 






^^#< 



:iA:::::/::/:/::: 1 i 

/: :\\ -.-.■/:; -J :f. ■/!.::>: 



V / ! A 






±i 



V: I 



■/■■■/ :- y- 



.-'• >-^ V 






<■ - 



sS± 



:/:::: 



:y ::: 



Figure 10. Example 2: The computed optical flow using the method of Horn & Schunck. 
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Figure 11. Example 2: The computed optical flow using the method of Cornelius k Kanade. 



21 
sped up by ignoring the offset term. We can instead use the updating equation 

to compute the optical flow and the multiplier fields. Here, we have defined 

u 



X a u - E x E t \ 
f'= ( v I, g'(f)= I X e v-E y E t 

rn t ) \\ m fnt + EEtJ 



and 



1 f E y X m + E X 8 + X m X s —E x E y X m E X EX 

(A')" 1 = -7 -E x E y \ m E 2 x X m + E 2 X 8 + X m X e E y E\ 



E X EX 3 EyEX a E 2 X X 8 + E 2 X„ + X* 



where 



a' = X m X 2 8 + E 2 X 2 + (E 2 X + E 2 y )X m X 8 , 

If the offset term is not negligible, then the estimates obtained from the above vector 
equation may be used as initial conditions in the original updating equations. 

7 Summary and Extensions 

Much of the existing methods for computing the local optical flow depend on two kinds of 
constraint: the flow field smoothness constraint and the brightness constancy constraint. 
The brightness constancy constraint permits one to match image brightness values across 
images. This constraint is sometimes very restrictive. 

We have proposed a new formulation by replacing the brightness constancy constraint 
with a more general constraint, which permits a linear transformation between image 
brightness values. The transformation parameters are allowed to vary slowly in space, 
so that inexact matching is allowed. We have formulated the problem of computing the 
optical flow as a minimization of a quadratic cost functional. Using variational methods, 
we have shown that the problem reduces to solving Laplacian equations for the two 
components of the optical flow field and the two transformation fields. We have described 
the implementation on a highly parallel computer, and presented sample results. 

One of the disadvantages associated with the use of the smoothness constraint is 
the degree of smoothness imposed on the unknown velocity and transformation fields. 
In fact, the algorithm developed here tends to smooth over discontinuities both in the 
optical flow and transformation fields. One way to overcome this shortcoming is to use 
"line processes" (Marroquin [1984]). Simply put, the idea is to incorporate in the penalty 
function the cost of introducing a discontinuity in the optical flow or transformation fields 
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instead of interpolating smoothly between two neighboring points when the gradient 
becomes large. A line process is a boolean variable; it takes on the value 1 where there 
is a discontinuity and 0, otherwise. 

There is, however, a drawback associated with the introduction of line processes in 
the minimization scheme, namely, that the cost functional becomes non-convex. This 
generally calls for inventing sophisticated optimization methods that can be computa- 
tionally exhaustive; for example, a simulated annealing scheme (Marroquin [1984]), or 
an algorithm based on neuronal network models (Hopfield &; Tank [1985], Koch et al. 
[1986]). 

Alternatively, we can employ the graduated non-convexity algorithm based on the 
concept of weak continuity constraints (Blake & Zisserman[l986]). Here, the process of 
minimizing a non-convex cost functional is replaced by minimizing a sequence of cost 
functions, the first of which is a convex approximation to the true cost functional and the 
last one is the true non-convex cost functional. Needless to say, we have yet to implement 
any of these schemes. 
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